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Abstract 

We investigate how Murray’s law is affected by fluids whose viscosity depends monotonously 
on shear rates, such as blood. Our study shows that Murray’s original law also applies to 
such fluids, as long as they did not undergo phase separation effect. When Fahraeus-like phase 
separation effects occur, we derive an extended version of Murray’s law. Finally, we study 
how these new laws affect the optimal geometries of fractal trees to mimic an idealized arterial 
network. Our analyses are based on Quemada’s fluid model, but our approach is very general 
and apply to most fluids with shear dependent rheology. 


1 Introduction 


Murray (1926) proposed the first law for the optimal design of blood vessels, based on a trade¬ 


off between the power needed to make blood circulate in the vessel and the metabolic power 
needed to maintain blood. His work was based on the original researches developed by [Hess] 


(1914). The law is formulated using Poiseuille’s regime in cylindrical vessels, thus it accounts for 


viscous effects only and neglects any perturbations due to fluid divergence or convergence at the 
bifurcations points. Blood flow rate is assumed constant to mimic for the fact that the vessels 
have to feed a downstream organ whose needs are independent on the vessels geometry. The 
optimal configuration corresponds to the blood flow rate being proportional to the cube power of 
the radius of the blood vessel. The well known corollary for a bifurcation states that the cube 
radius of the parent vessel equals the sum of the cube radii of the daughter vessels. Murray’s 
law is independent on the amount of blood flow rate, thus it does not depend on the functioning 
regime of the downstream organ, at least in the limit of its hypotheses. In the seventies, [Zamir 
(1976) extended this law to account for bifurcation angles and he expressed Murray’s law in term 
of wall shear stress being constant independently on the vessel size (Zamir 1977). While blood 


arterial mid-level circulation meets conditions for Murray’s law, this is not the case for the larger 


vessels, where inertia and/or turbulence occur ( 

Uylings 

L977 

), and for microcirculation, where 

wall shear stress is decreasing with the sizes of the vessels 

(Sherman 

1981). Different hypotheses 


were developed to explain the shifts to Murray’s law observed with microcirculation, for example 
Taber (19981 proposed to add an energy cost related to smooth muscles. Later 


(2005) used semi-empirical laws from Pries et al. (1994) and showed that wall shear stress behavior 


Alarcon et al. 


in microcirculation can be explained by the Fahraeus-Lindvquist effect. The Fahraeus-Lindvquist 
effect is a phase separation effect occurring in small blood vessels that makes blood viscosity 


become a non monotonous function of vessel radius (Fahraeus 1929 Pries et al. 1995 Fung 1997 


Pries and Secomb 2008). Similarly, studies were made to understand why Murray’s law does not 


apply in large blood vessels, where inertia and turbulence play an important role. These studies 
were based on modified formulation of the power associated to fluid circulation, shifting the radius 
exponent from 3 in laminar case down to 2.33 for fully turbulent flow (Uylings 1977). Finally, 
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generalization of Murray’s law has also been developed to account for the non-local properties of 
tree structures and extends the predictions using empirical exponents based on the fractal nature 
of the biological networks (Zhou et al. 1999 Kassab 2007). These exponents aggregate rich 


information from the physiological configuration, amongst which the non linear behavior of blood 
viscosity. These last studies should however be considered with care, since the empirical data 
injected in the model is actually dependent on the variables that are optimized, which might bring 
a bias in the predictions. 

Another step to fully extend Murray’s law to blood circulation is now to add the dependence 
of vessel blood viscosity on flow amplitude. Indeed, blood is a shear-thinning fluid and its local 
viscosity depends on the local shear rate in the vessel. Local shear rate is a function of both blood 
flow amplitude in the vessel and vessel radius. Thus, in addition to red blood cells volume fraction 
and vessel radius, equivalent viscosity in blood microcirculation is expected to be also dependent 
on the shear rates in the vessel. 

First we show that if Fahraeus effect does not occur, the fluid equivalent viscosity in a vessel 
is dependent on mean shear rate and hematocrit in the branch only. If Fahraeus effect occurs, 
then we show that equivalent viscosity becomes also dependent on vessel radius. Then we apply 
Murray’s optimal design to both cases with and without Fahraeus effect and derive corresponding 
laws for optimal configuration of a vessel and of a bifurcation. Finally, we study how these new 
laws affect optimal geometries of fractal tree structures, using inherent properties for mean shear 
rates variation inside a fractal tree. 


2 Mathematical model 


Since Murray’s law was formulated in the frame of the cardiovascular system, we chose to work 
with a model that mimics blood rheology. However, blood exhibits a complex thixotropic behavior 
with a yield stress to overcome for blood to flow (Merrill 1969[ Bureau et al. |1980 Apostolidis and 


Beris, 2014). Moreover, blood rheology is not only affected by red blood cells concentration but 


also by its inner composition in proteins, such as fibrinogen concentration (Merrill et al. 1969). 


We did not want to account for such refined behaviors of blood rheology since our study focused 
on the interaction between Murray’s optimization process and a shear rate-dependent rheology. 
Thus we chose to work with Quemada’s fluid model which tightly fits our needs and were initially 


proposed for modeling blood rheology. Those fluids are well documented (Quemada 1984) and 


have been intensively studied and validated as good approximations for blood modeling (Cokelet 
1987 Neofytou 2004 Marcinkowska-Gapinska et al. 2007 Sriram et al. 2014). The viscosity in 


Quemada’s fluid model depends on the local shear rate 7 and on the local red blood cells volume 
fraction H. Complete details about Quemada’s model are given in section I of supplementary 
materials. Notice that other models, such as Casson’s model could also have been used in this 
study ( |Casson , |1959| . 

In Quemada’s model, the dependence of fluid viscosity on shear rates divides into three parts, as 
shown on figure [T] (blood case): a plateau of high viscosity for low shear rates (blood: < 10 -3 s _1 ); 
one sharp decrease of viscosity for ’’medium-ranged” shear rates (blood: from 10~ 3 s _1 to 1 s _1 ); 
a plateau of low viscosity for high shear rates (blood: > 1 s _1 ). 

We assume the vessels to be cylindric and the fluid to be at low regime, axi-symmetric and fully 
developed in all the vessels. Blood flow rate through a vessel is assumed constant in order to mimic 
a downstream organ which is fed by the vessel and whose needs are independent on the delivering 
structure. Then combining Quemada’s formula for viscosity with fluid mechanics equations in a 
vessel, we use the equivalent viscosity fi eq of our fluid in the vessel. It is defined as the viscosity a 
Newtonian fluid would need in order to dissipate the same amount of viscous energy than our fluid 
inside that vessel. Interestingly when Fahraeus effect can be neglected, the equivalent viscosity 
is driven only by the mean shear rate in the branch ( 7 ) and by the red blood cells volumetric 


fraction Hp: ^ eq (( 7 ) ,-ffo) = > with G a smooth function. Explanations about why and 

how they are such related are detailed in section II-A of supplementary materials. Fahraeus effect 
occurs mainly for vessels whose diameters are smaller than 300 \wn, see for example [Pries and 
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Figure 1: Quemada’s model of viscosity (Quemada 1984 Cokelet, 1987). Black curve: viscosity 
dependence on the shear rate (blood case, H = 0.45). Red dashed line: optimal mean shear rate 
predicted by the model for blood large circulation ( 7 )* noF ~ 12.5 s -1 . 


or section III of supplementary materials. For such vessels, the independence of 
viscosity on the vessel radius r is lost. As a consequence, the equivalent viscosity 
can be expressed using a smooth function K: n eq ({ 7 ) ,r,H D ) = , see details in section 

II-B of supplementary materials. 


Secornb (2008 


the equivalent 


3 Extending Murray’s law to Quemada’s fluid model 

In this section, we extend Murray’s law to non-Newtonian fluids that can be modeled with 
Quemada’s fluid model. Thus we take into account the dependence of viscosity on shear rates. 
We start with the assumption of no phase separation effects in the fluid. Next, we study the case 
where Fahraeus effect occurs in the fluid. 


3.1 Without phase separation effects 


Now we will show how these principles allow to extend Murray’s law to any fluid that can be mod¬ 
eled with Quemada’s model -and consequently to blood- as long as Fahraeus effects are negligible. 
Let us consider a fluid that can be modeled with Quemada’s model, and let us assume its flow 
rate F is going through a vessel with radius r and length l. The dissipated power W defined by 
Murray ( |Hess 1914 Murray 1926 Alarcon et al. 2005) divides into two parts: W — Wh + Wm, 
where Wh is the power dissipated by the flow, and Wm the energy consumption rate of the fluid 
(in the case of blood, this is a metabolic consumption rate). With the above results and denoting 
ab the energy consumption rate per unit volume of fluid, we have 


W H = 


8FV 9 «7))Z 


and Wm = otb'Krl 


(i) 


The design principle proposed by Murray is to search for a minimum of W relatively to the radius 
of the vessel, thus solving ^ = 0. In our case, expanding this last equality leads to a non linear 
equation which depends only on ( 7 ), see section III of supplementary materials. As a consequence, 
the optimal configuration is reached when the mean shear rate in the vessel solves the preceding 
equation, independently on the fluid flow rate F or the vessel sizes r or l. Interestingly, we 
observed that the optimal shear rates are as small as possible, to get low viscous dissipation, but 
still high enough to keep viscosity in its lowest values. In the case of blood, the optimal shear rate 
is (7 )* no F ~ 12.5 a- 1 , it can be easily located on figure [l] it stands on the left of the low viscosity 
plateau (red dashed line). 

Let us apply this result to a bifurcation where the flow rate in the parent branch is F p and the 
flow rates in the daughter branches are F\ and F^- Their respective radii are denoted r p , n and r 2 - 
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If (j)n oF i s the solution of ^ = 0, then ( 7 )* noF is the optimal mean shear rate in the three vessels 
and (7 )* no F = ^3 = = By flow conservation, we have an additional equation which is 

Fp = F-| + F'2- Combining these equations finally leads to Murray’s original law: r 8 = rf+r 3 . 

The method and these result are most general, they apply to any fluid whose viscosity is 
monotonously driven by the shear rate in the vessel. 


3.2 With Fahraeus effect. 


The Fahraeus effect in blood is a phase separation effect due to the biphasic characteristics of 
blood: red blood cells tend to migrate toward the centre of the vessels and blood near the vessel 
wall is depleted in red blood cells. Fahraeus effect contributes to blood Fahraeus-Lindvquist effect 


(Fahraeus 1929 

Fung 

1991 

7 ). It becomes non 

negligible for vessels with diameters smaller than 

300 /im, see for example ( 

3 ries and Secornb 

2008 

). To estimate the role of Fahraeus effects 


on Murray’s law in such vessels, we approximated this effect by assuming that a red-blood-cell 
depleted layer stands near the wall of the vessels, see section II-B of supplementary materials. The 
thickness of this depleted layer, and consequently the whole fluid dynamics into that vessel, depend 
on the branch radius (Pries and Secornb, 2008). As a consequence, the equivalent viscosity in a 
branch where Fahraeus effect occurs is not anymore dependent on the mean shear rate only, it also 
becomes dependent on the radius of the branch: n eq (( 7 ), r) = , with K a smooth function. 

In this case the power dissipation in the fluid used for deducing Murray’s law is Wh = 8F . 

The energy consumption rate of the fluid Wm remains Wm = ctbFr 2 l. As before, the radius that 
minimizes the total work W = Wh + Wm solves = 0. Expanding this last equation shows 
that the minimum power is reached on a curve r —► (q) F (r) that depends on the red blood cells 
volumetric fraction FLp>. As an example, we plotted the curve computed in the case of blood on 
figure © When r is large enough, say larger than about 300 /jm ( |Pries and Secomb| |2008|), then 
the dependence on r is lost and ( 7 )^, (r) = ( 7 )* oF - 



Figure 2: In black, the optimal mean shear rate in a branch as a function of the radius of the 
branch with Fahraeus effect (function r —> ( 7 )* noF (r) in text); in red dotted line, the optimal mean 
shear rate without Fahraeus effect. Blood case: = 0.45 and ab = 77.8 J.m~ 3 .s~ l (Taber 


1998 Alarcon et al. 2005). 


Let us now consider a bifurcation where the mean shear rate in the parent branch is ( 7 ) and 
the mean shear rates in the two daughter branches are ( 7 ^ and ( 7 ) 2 - Their respective radii are 
denoted r p , r\ and r?- Using flow conservation through the bifurcation, we can deduce that the 
optimal radii occur for the following equality: 

(7 )*f Op) 4 = (7 }*f 00 r i + (7 )* F O 2 ) r\ (2) 

As expected, when Fahraeus effect vanishes, the previous equation simplifies into the original 
Murray’s law r 8 = rf + r 8 . 
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4 Extended Murray’s law in fractal trees 


We are now interested in a tree structure built as a cascade of cylinders. The branches divides 
regularly into n smaller identical branches. The number of divisions between a vessel and the 
root vessel of the tree defines its generation index. The tree root stands at generation 0 and the 
tree limbs end at generation N. The size of the branches are defined thanks to the homothety 
ratio or scaling factor h (Mauroy et al. 2004 Mauroy and Bokov, 2010) that corresponds to the 
relative change in vessels diameter and lengths in a bifurcation, i.e. if ri is the radius of vessels 
of generation *, then the radii of vessels of generation i + 1 are r i+1 = h x 77 . It has been shown 
that bifurcating points in tree structures affects in a complex way the fluid dynamics, even for 
the Newtonian case, see for example Pedley et al. (1970). However, for the sake of keeping the 
model tractable, we will stick to Murray’s initial hypothesis and assume that fluid dynamics is not 
perturbed at the bifurcations and that velocities remain fully developed all along the branches. 
Blood flow rate through the tree is assumed constant to mimic the needs of a downstream organ 
that are independent on the network geometry. 



Figure 3: Tree network structure with n = 2 and N — 6. The tree is dichotomous and the vessels 
size decreases at each generation: their diameters and lengths are multiplied by the homothety 
factor h < 1 after each bifurcation. 


4.1 Mean shear rates 

Denoting F, the blood flow in a branch of generation i and Si = irrf the surface of its circular cross 
section, then the mean shear rate in that branch is ( 7 ,;) = F, ^ S ' = 4-y. Since the tree branches 
divide into n smaller identical branches, the total blood flow in a branch of generation i is n times 
the total blood flow in a branch of the next generation i + 1. Then the mean shear rate in a branch 
of generation i + 1 is ( 7 ^+ 1 ) = = ( 7 ») T J 4 » • Thus, the regularity of the structure 

(fractal) induces also a scaling law on the mean shear rates in the tree. This scaling law depends 
only on the homothety ratio h. Depending on the position of the factor — 7-3 relatively to 1, the 
mean shear rate has different behaviors: 

l/o 

1. If h > (1) , then the mean shear rate decreases along the generations, consequently blood 

viscosity tends to increase along the generations; 

■l/o 

2. If h = (4) , then the mean shear rate remains constant along the generations and so for 

blood viscosity; 

1 /q 

3. If h < (4) , then the mean shear rate increases along the generations, consequently blood 

viscosity tends to decrease along the generations. 
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Figure 4: Mean viscosity variation in a dichotomous tree (n = 2) with TV = 10 generations for 
three different values of the homothety reduction factor h. In the case plotted Fahraeus effect is 
negligible and the shear rate in the root branch of the tree is 640 s _1 . 


Blood viscosity depends on shear rate in a non linear way, with a plateau of high viscosities at 
low shear rates and a plateau of lower viscosities for high shear rates. Between the two plateaus, 
for medium shear rates, blood viscosity varies steeply, see figure [l] When the shear rate decreases 
or increases along the generations of the tree, then viscosity will vary more strongly if shear rate 
variation makes the viscosity go through the steep part. Consequently, the amplitude of viscosity 
variation depends on the initial mean shear rate ( 70 ) in the root branch of the tree: a high shear 
rate in the root branch can induce notable viscosity variations throughout the tree only if shear 
rate is decreasing enough along generations; similarly, a low shear rate in the root branch can 
induce notable viscosity variations throughout the tree only if the shear rate increases enough 
along the generations. In any other case, the viscosity is stalled, either in one of the two plateaus 

1 /3 

or because mean shear rate does not vary much if h ~ (4) 


4.2 Without phase separation effects 


In the absence of Fahraeus effect, Murray’s law extended to Quemada’s model states that the 
mean shear rates in the tree branches should be all equal to the optimal mean shear rate ( 7 )*. 
Since the mean shear rate is following a scaling law in the tree, the only way for the tree to 
minimize the dissipative power throughout the tree is for the scaling parameter —p- to be equal 

to 1, i.e. Ti* = ( 7) 3 • This optimal configuration corresponds to mean shear rates and equivalent 
viscosities being constant throughout the whole tree. This result is fully compatible with blood 

2008). For blood, it has been estimated that otb = 

with these numbers, 


arterial macrocirculation (Pries and Secomb 
77.8 J.m~ 3 


\s 1 and Hd = 0.45 (Taber 


optimal mean shear rate we predict is ( 7 ) 


1998 


Alarcon et 
=T 


al. 


noF 


12.5 


see 


j 2005 ) 
figure [I 


the 


Interestingly, except in very large vessels and in microcirculation, blood arterial network ge- 

with 


2008) 


ometry exhibits a constant mean s hear rate along its generatio ns (Pries and Secomb 

a homothety ratio of about (^) 3 (Rossitti and Lofgren 1993), as predicted by our model. This 
aises the question on how physiology can reach this configuration. A scenario proposed by Zamir 


(1976) is based on the today well-known fact that endothelial cells standing on the arterial walls 
are able to respond to shear stress stimuli (Pries and Secomb 2005). If this response is made 


accordingly to a shear stress threshold common to all cells, increasing the vessel diameter if the 
shear stress is over a threshold (by cells division) and decreasing the vessel diameter if the shear 
stress is below the same threshold (either by apoptosis or migration), then the resulting tree will 

exhibit a homothety ratio that is (^) 3 • This is a direct consequence of the scaling law on mean 
shear rates described in the previous section. 
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4.3 With Fahraeus effect. 


The case of arterial microcirculation is more complex, because it includes phase separation effects. 
Many experiments have shown the influence of Fahraeus effect on equivalent viscosity, as reported 


Pries et al. (19921. Fully detailed analysis and curve fittings are available in Pries et al. (19921 


and Pries and Secomb (2008) about the relative viscosity, which is the ratio between the equivalent 
viscosity and the embedding fluid viscosity (plasma in the case of blood). Our model remains in 
the range of the experimental data and gets a similar trend as that of the data in these studies, 
as shown on figure [5j Notice that by altering the parameters of our model, we would be able to 
better agree with the fit given in Pries et al. (1992), but this is out of the scope of this paper. 



Figure 5: Relative viscosity versus tube radius for a constant pressure drop in the tube (high 
shear rates). The dashed red line represent the mean fit from experiments performed in Pries 


et al. (1992). The black continuous line represents the behavior of our model. Notice that we are 


interested only in having a similar trend and not to exactly fit the red dashed curve. 


The optimal configuration of each division in the tree meets the conditions for the extended 
Murray’s law equation ©> with daughter branches all identical. In addition, the optimal con¬ 
figuration also verifies the scaling law for the branches radii, i.e. ?y + i = /i?y. These equations 
make the scaling factor h generation dependent, thus scaling factors need now to be indexed 
with their generation index i. The optimal hi in term of Murray’s design verifies the equation: 
nh l* (h it * r i,*) = (7)f ( r h*)- 




Figure 6: Example for blood circulation in a dichotomous tree (n = 2), with 9 generations, a root 
radius of 800 /im and a rati o of length over radius for the branches equal to 10; Blood properties 
are from (Taber 1998; Alarcon et al. 2005): Hjj = 0.45 and ab = 77.8 J.m -3 .s^ 1 . (a): Optimal 
scaling factor with Fahraeus effect (black) and without Fahraeus effect (red), (b): Mean shear 
stress versus pressure within the tree. Deepest generations are on the left and upper generations 
are on the right. 
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We computed numerically the optimal scaling factors in the case of blood for a dichotomical 
tree (n = 2) with nine generations, see figure [6j4. The root branch radius is 800 /im, thus Fahraeus 
effect is negligible in the first generations. In the optimal configuration, the three first generations 

are bifurcating with a scaling factor of (1 ) 3 ~ 0.7937, then Fahraeus effect becomes non negligible 
and it affects the optimal scaling factors of the next generations. The function r —» ( 7 )^ (r) is a 

decreasing function, thus /ij ; * < (^) 3 = h*. Thus the optimal configuration has tighter branches 
where Fahraeus effect occurs than the optimal configuration without Fahraeus effect. Indeed, 
Fahraeus effect improves the lubrication of the core fluid in the tube because the red blood cell 
depleted layer near the wall is less viscous that the core layer. Consequently, equivalent viscosities 
decrease and mean shear rates increase along the generations as soon as Fahraeus effect appears. 

The mean shear stress in a branch where Fahraeus effect occurs is 07 * = ^K((-y) F (r^*), r^*). 
Since K is an increasing function of the mean shear rate, the mean shear stress increases along the 
generations when Fahraeus effect occurs (see section II-B of supplementary materials for details). 
Because pressure is decreasing along the generation, the mean shear stress is a decreasing function 
of pressure within the tree, see the example of blood on fig ure [ 6 )>. This r esult is in agreement 


with the observation of arterial blood microcirculation (Pries and Secornb 2008), and confirms 


the conclusions of Alarcon et al. (2005) that Fahraeus effects might be the core phenomena driving 
the decrease of blood equivalent viscosity and wall shear stresses in the microcirculation. 


5 Conclusion 

In this study, we extended Murray’s law to fluids that can be modeled with Quemada’s model. 
In a cylindrical tube, we showed that the mean shear rate drives the behavior of the equivalent 
viscosity of the tube and thus of the power dissipated in the flow. The consequence is that Murray’s 
optimization principle in a cylindrical tube brings an universal optimal mean shear rate that does 
not dependent on the flow amplitude or on the size of the tube. In the case of Quemada’s law 
mimicking blood behavior, we found that this optimal mean shear rate is ( 7 )* noF ~ 12.5 s -1 . 
This value is in accordance with arterial macrocirculation data where the mean shear rates in 
the successive vessels remain almost constant at about 10 s -1 . Applied to a fractal tree whose 
branches divide into n identical smaller branches, we found that the optimal tree scaling factor 

JL 

remains the same as for a Newtonian fluid, i.e. (^) 3 . 

However, if the vessel diameter is small enough, Fahraeus phase separation effect occurs. A red 
blood cells depleted layer appears near the wall of the vessel, thus making the equivalent viscosity 
dependent on the mean shear rate and on the vessel radius r. We derived a new optimal config¬ 
uration from Murray’s optimisation principle, which is expressed through a decreasing function 
that relates the optimal mean shear rate in a vessel with its radius: r —> ( 7 )* F (r). For large r, this 
function maps the optimal mean shear rate without Fahraeus effect. The function behavior relates 
to the lubrication phenomenon induced by Fahraeus effect. We also derived the optimal configu¬ 
ration of a fractal tree when Fahraeus effect occurs. We showed that the optimal scaling factors 
become dependent on the generation index, on the size of the root of the tree and on the function 
r —> ( 7 ) F (r), i.e. nh %* ( / y) F = ( 7 )^ ( 77 *). The optimal configuration of the fractal tree 

allows tighter branches when Fahraeus effect occurs and induces a drop in mean shear stresses 
along the tree, which has been reported in the literature for blood arterial microcirculation. 

Our study is based on the sole fact that viscosity is a monotonous function of the fluid shear 
rate. Consequently, our results hold not only for fluid that can be modeled with Quemada’s 
model, but also for any fluid for which this monotonous condition is true. This make our study 
most general and apply for example to shear-thinning or shear-thickening fluids. Concerning 
blood, interestingly, the optimal configuration fits that of large circulation. Large circulation is 
thus in a state where non-Newtonian effects remain small, just on the edge of the steep part of 
blood rheogram, as shown on figure [lj It is important to notice that this property is a result of 
the optimization process and not an hypothesis. This result can only be obtained by including 
the non Newtonian behavior of blood in the optimization process. This also highlights why, 
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although it was based on Newtonian hypothesis, the original Murray’s law fitted nevertheless blood 
large circulation. Of course, many other phenomena might play a role in the complex behavior 
of arterial blood circulation. Hence, our model does not account for the complex regulation 
processes that exist in blood network. These regulation processes are likely to affect the optimal 
configuration and enhance the robustness of the system. In large circulation, the role of inertia, 
turbulence or oscillating flow caused by the heart beat may influence Murray’s optimal design. In 
the microcirculation, the roles of other phase separation effects have to be taken into account to 
improve the validity of the predictions, for example the role of the glycocalyx molecules standing 
on the vessels walls (Pries et al. 2000). Also, because the power dissipation in Murray’s optimal 
design is not symmetric relatively to its optimal value, biological noise affecting blood vessel sizes 
may also influence the optimal configuration (Mauroy and Bokov 2010 Vercken et al. 2012). The 
inclusion of such a noise in Murray’s optimal design may bring interesting insights on blood vessels 
size distribution in the arterial network. Finally, blood network is regulated and this regulation 
probably affects the optimal configuration. 
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Appendix 


A Quemada’s law 


The local viscosity /i is a function of both local shear rate 7 ( 5 ) and local hematocrit H(s) as stated 
by Quemada’s law: 


Mi,#) = Fp 


h y 2 ( i +k(H) \ 

H^H)) \ X (H) + k(H)J 


1 - 
H 0 and 77, 


H 

HUH) 


with k{H) = (t c (#)lil) 2 > X(H) = (l - h)(h)) / ( 

/i p = 1.6 10 -3 Pa.s. The quantities t c (seconds), L 
volumetric fraction H. We use the fits from (Cokeletl 1987): 


)\ is the viscosity of plasma and 
are functions of the red blood cells 


^ _ g 6.1508—27.293if+25.60ii 2 —3.697ii 3 


H 00 (H) = 2e- 


1.3435+2.80344-2.7im 2 +0.6479.F/ 3 


Hq(H) = 2 e - 3 - 8740 + 10 - 410 #- 13 - 80 # 2 + 6 - 738ff3 


The formulas for f c , and Hq come from data fitting, they are valid for red blood cells 
volumetric fraction up to 0.7 (Comolet 1984) and for moderate shear rates, i.e. larger than 
about 0.1 s -1 . Other fitting methods using second order polynomials can be used, see for example 


(Sriram et al. 2014). Using such alternative fit induces very few differences on the results predicted 


in this work. 


B Computation of the effective viscosity 

We assume the vessels to be cylindric and the fluid to be axi-symmetric and fully developed in all 
the vessels. Following Quemada’s law, the viscosity /i (7 (s),H( s )) of the fluid at a radial position 
s of a vessel depends on the local shear rate 7 (s) and on the local red blood cell concentration 
H(s). The pressure drop per unit length is denoted C, it is assumed independent on the position. 
In this frame, the fluid dynamics reduce to an equation on the shear rate 7 (s) at radial position 
s e [0, r]: n('y,H)j = 

Once the expression of the viscosity /4 is known, this equation shows that the shear rate at 
radial position s, 7 , is a function of Cs and H only. Using Quemada’s law, it is possible to reach 
an analytical expression for the shear rate 7 , and consequently for the flow rate F in the branch 
as a function of the radius of the branch r, of the pressure drop per unit length C and of the 
discharge hematocrit Hjj, see appendix [B] for the details. Note that for branches whose radius is 
sufficiently large, say larger than 300/im, Fahraeus effect is negligible and the blood hematocrit H 
in the branch is homogeneous and equal to the discharge hematocrit Hn>, see figure [ 8 ] From these 
relations, and the expression of the hydrodynamical resistance of a cylindrical tube of length L , 
R = CL/F = ( 8 /<L)/( 7 rr 4 ), we can compute the effective viscosity )i eq of the fluid in a branch 
depending on the quantity Cr in the branch and the red blood cells volumetric fraction FLd in the 
branch: /r eg (r, C, H D ) = sF(r,c°H D ) = 8 ( 7 )(j U,Hd) ' 

B.l Without Fahraeus effect. 

In this section, we consider a cylindrical branch with radius r and length L. Blood flow in the 
branch is assumed to be fully developed and axi-symmetrical. The quantity C is the pressure 
drop per unit length in a vessel, it is assumed to be a constant in the whole vessel. We do not 
take Fahraeus effect into account, ths red blood cells volume fraction H is assumed constant in the 
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whole vessel. The shear rate at radial position s is then denoted j(s, C , H) and the fluid dynamics 
in the branch with these hypotheses reduces to and is the solution of the following equation: 

fa 

A* ( 7 («, C, H), H(s)) 7 (a, C,H) = — ( 3 ) 

The first consequence of this equality is that the shear rate y(s, C, H) can actually be seen as 
a function of Cs and H only : 7 (Cs,H). The quantity Cs is a constraint; denoting x the axis 
of the branch, 2 irCsdx represents the infinitesimal increase of the pressure force felt by a particle 
of fluid whose radial position in the branch is s and which moves a distance dx forward into the 
branch. 

The second consequence of equation (J3J) is that we can compute an analytical solution for the 
local shear rate 7 (Cs,H): 


( 


i(Cs,H) = 


yAyi (Cs,H) 


y/i 3 (H) ^(V7 i( C s , h ) ~ +^'ii{Cs,H)^ 2 {H) 


V 7 

(4) 

with 71 (Cs,H) = ^ (l - > 72(-ff) = and 73(F) = Then the flow F in a 

branch is: 


F = 


r s 2 

- 27 r / j(Cs, H)—ds 

Jo 2 


( 5 ) 


Consequently, the mean shear rate in the branch (7) is a function of Cr and H , as shown by 
the following equalities: 


( 7 ) = 


F 

7 T r 3 


fo 7 (Cs,H)s 2 ds _ f 0 Cr j(y,H)^±dy 


/ 0 Cr 7 (y, H)y 2 dy 
(Cr ) 3 


( 7 ) (Cr,H) 


( 6 ) 

Moreover, the mean shear rate (7) is a strictly increasing function of the pressure drop per 
unit length (and as a matter of fact of the radius of the branch too, since it depends on Cr), thus 
there exists a function G such that Cr = G(( 7) ,F) that relates the mean shear rate (7) to Cr 
and H. No analytical expression for G is available but it can be easily computed numerically by 
inverting equation ([6]). 

From these results, we can get an analytical formula of the effective viscosity of the fluid in a 
branch as a function of Cr and H: 


y eq (Cr,H) 


7rr 4 C 
8 F 


or using (7) as the variable: 


1 Cr 
(7) (Cr,H)T 


Cr A 
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7 (Cs,H) —ds 


-1 


( 7 ) 


deq (( 7 ) ,H) 


l g((7),g) 

8 (7) 


( 8 ) 


To summarize, when Fahraeus effect can be neglected, the quantity Cr is a function of the 
mean shear rate in the branch (7) and of the red blood cells volumetric fraction Hjj, they are 
related through the function G: Cr = G((-j) ,H D ). Thus, when Fahraeus effect does not occur, 
then the effective viscosity is driven by the mean shear rate in the branch (7) and red blood cells 
volumetric fraction He>: y eq (( 7) ,Fd) = • 
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B.2 With Fahraeus effect. 


In this section, we consider a cylindrical branch with radius r and length l. Blood flow in the 
branch is assumed to be fully developed and axi-symmetrical. The quantity C is the pressure drop 
per unit length in a vessel, it is assumed to be a constant in the whole vessel. Blood hematocrit 
is not anymore assumed constant in the whole vessel since we account for Fahraeus effect in this 
section. 

Red blood cells distribution in a vessel. 

We assume that red blood cells stand only in the core of the vessel where the hematocrit is 
constant and equal to H (core hematocrit). The radius of the core is £ < r. Plasma flows in the 


considered, that are larger than 50 microns. Plasma and red blood cells flow in the core (their 
flow is F core ), see figure [7j 


layer near the wall (its flow is Fi ayer ). We neglect the role of the glycocalyx layer (Pries et al. 
2000 ) since its thickness, of about 1 micron, is small relatively to the diameters of the vessels 


RBC hematocrit: 0 

RBC hematocrit: H 


layer 


core 


Figure 7: Scheme of the blood flow inside an axi-symmetric vessel. 

The red blood cells depletion near the vessel wall induces the Fahraeus effect that corresponds 
to a decrease of the mean hematocrit in the vessel relatively to the discharge hematocrit Ho in 
large vessels. Since the blood that circulates into the vessel comes from larger vessels, by flow 


conservation (Fung 1997 Pries and Secomb 2008) we can relate discharge hematocrit Flo to core 
hematocrit FI 


H d = 


H{r, C, H d ) F core (r,C,H D ) 


(9) 


Fcore (r? C. FIo ) T Fi a y er (/'. (F • Ho ) 

The dependence of the mean red blood cells volumetric fraction in a vessel, called tube hemat¬ 
ocrit Ht , with tube radius and discharge hematocrit is well documented (Pries and Secomb 2008): 
if the radius of the tube r is measured in microns, then 


H T (r, Ho) =H d x ( H d + (1 - H D ) x (1 + 1.7 e -°- 415x2r - 0.6e-° ollx2r ) 


( 10 ) 


Finally, tube hematocrit is related to core hematocrit by Ht{t, Ho)^r 2 = H(r, C, Ho)n£(r, C, Ho) 2 


and 


H T (r, H d ) = H{r, C, H D ) ^ (11) 

We can thus define the following function, which represents how red blood cells volumetric 
fraction varies along the radius of the vessel: 

if 0 < s < £ (core), H(s) = H{r,C,H D ) = ^ r ^ Hn)2 H T (H D ,2r) 

if £ < s < r (wall layer), H(s) = 0 

In order to be able to compute the values of H(r,C, Ho) and £{r,C, Ho), we need to 
this set of equations the equations of fluid dynamics. 

Fluid mechanics. 

The shear rate 7 at radial position s is the solution of the following equation: 


( 12 ) 
add to 
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r's 

/i(i(8,C,H(s)),H(s))'y(s,C,H(s)) = — 


(13) 


As in the case with no Fahraeus effect, this shows that shear rate at radial position s is a 
function of Cs and H{s). However, on the contrary of the case without Fahraeus effect, the 
dependence of 7 on s remains through the dependence of H(s) on s. From equation (13), we can 
compute an analytical solution for the local shear rate 7 {Cs, H(s)): 


if 0 < s < £ (core), 

_ ( \/7i(Cs,ff(s)) y/-y 3 (H(a)) , \j (V'Vi( C ' s - H ( s ))-V' t's^W)) ( ff (»)) 

< 7(0 s, JrL fsjj — I 2 2 ' 2 

if £ < s < r (wall layer), 7 (Cs,H(s)) = 

(M) 

with 71 (Cs, ff(s)) = jgj (l ~ huhU) ) ’72( h (s)) = £§^£ and 7 3 (ff(a)) = 

As a consequence of the dependence of the shear rate on s, the mean shear rate and the effective 
viscosity in the branch are both dependent on the radius of the branch r: 


( 7 ) (r,C,H») 


Fcore{^i C, H £)) T Flayer ^A C, H-d) 
7rr 4 


(15) 


IJ j eq(j‘> 1-0 Hd) 


Cr 

8(7) (A <?,#£>) 


(16) 


Once r and -Hd are fixed, the mean shear rate ( 7 ) in the branch is a strictly increasing func¬ 
tion of the pressure drop per unit length C, thus there exists a function K such that Cr = 
K{{ 7 ) The function K is computed numerically (see next section). Thus, effective vis¬ 

cosity can be reformulated into new variables: 


/“eg((7) i r,H D ) 


Kid) ,r,H D ) 
8 ( 7 ) 


(17) 


Solving the system. 

In this section, we first explain how the diameter of the red blood cells core £ and the volumetric 
fraction of red blood cells in the core, H , are computed once we know the radius of the branch 
r, the pressure drop per unit length in the branch C and the discharge hematocrit Hjj. Then we 
explain how we obtain the function k that relates C to the mean shear rate in the branch ( 7 ). 

The red blood cells core £ and the volumetric fraction of red blood cells in the core H are the 
solution of the following system that regroup the conservation equations in the branch and the 
equations from fluid dynamics (see previous sections): 


u _ u F core (r,C\Hn>)-\-Fi ayer {r,C\H £>) 

n ~ 110 F cor e(r,C,Hr>) 

c2 — H T {r,H p ) „2 
? ~ H 1 

F core (r,C,H D ) = 27r (v(£)^- - j(Cs, H(s))^-dsJ 

Flayer (A C,H D ) = -27T (jj j(Cs, H (s)) ^ds + u(£) ^ 

u (0 = ( r2 -£ 2 ) 

7 (Cs,H(s)) is given by equation Q 


(equation ®) 

(equations ( flO] ) and @) 


(18) 


The system (18) is normalized and solved using Matlab, thanks to the fsolve command (based on 
a Newton method). To compute the two integrals in equation (18), we use the quad command 
(adaptative Simpson quadrature). 
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The mean shear rate in the branch ( 7 ) ( Cr,Hz j) is then computed using equation (151. The 
function k that relates the quantity Cr to the mean shear rate ( 7 ) is then computed numerically 
by numerical inversion of the function ( 7 ) ( Cr,Hu ). 





Figure 8 : A: Fahraeus effect plays a role on branch with small radii only, typically smaller than 
300 \im. B: Effective viscosity /i eq depends on mean shear rate in branch unless Fahraeus effect 
occurs, case with mean shear rate of 2/n x 10 ~ 3 s -1 . C: Dependence of effective viscosity p, eq with 
mean shear rate (example with tube radius 200 fj,m where Fahraeus effect occurs and plays a role 
mainly at low mean shear rates (< 0.1 s” 1 )). 

To summarize, when Fahraeus effect occurs, the results computed without Fahraeus effects 
remain a good approximation if the vessel radius is larger than about 300 /rm. For smaller vessels, 
we lose the independence on the radius of the branch r, and the effective viscosity can be expressed 
using a function K relating Cr to ( 7 ): /r eq {(i ) ,r,Hu) = = K ^^Q } Hd ' > . In general, Hu is 

constant and it will dropped off the equations for clarity reasons. All our results would nevertheless 
depends on the value of H D . 

Finally, the mean shear stress in a branch is ay* = M eq{{i)* F (r i ,*), ry*) ( 7 )^ (ry*) = |IF(( 7 )^ (ry*), ?y*)- 


C Optimisation of Murray’s cost. 


The numerical example we gave in this work for optimal mean shear rate use Quemada’s fluid 
and the set of parameters corresponding to blood (see append ix |A|. They were computed with 
Hu = 0.45 and ab = 77.8 J.?n _ 3 .s _1 (Taber 1098 Alarcon et abj 2005). The numerical method 
depends on whether Fahraeus effects occur or not. The first step is always to build the function 
relating the effective viscosity to mean shear rate, hematocrit and, in the case where Fahraeus 
effect occurs, with vessel radius, see appendices [B| 

= 0 , i.e. 


1 dW 


When Fahraeus effect does not occur, then we solved directly the equation ^yyyyr 

&l*Leq 


nrl dr ^ 


4 m(( 7 )) + 3 ( 7 ) 


5 ( 7 ) 


«7» + 2a b = 0 
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This formula shows that the minimum of W is reached for a given value of (j)* loF - The derivative 
is computed numerically using a first order method. 

When Fahraeus effect occurs, W is given by 


W 

nl 


8 F 2 ^ eq ({j) ,r) 

'77-2 r*4 


+ a b r 2 


and its derivative of W relatively to r expresses as 


7 = ( 4/X + 3 ^ ^^0 + 2<Xb = 0 


This formula shows that the minimum power is reached on a curve r —> (j)* F (r). In the case 
with Fahraeus effect, we did not work with but instead we minimized directly the cost W = 
Wh + Wm using a gradient method. A value for the fluid flow rate F is chosen as stated in 
Murray’s optimal design, then we compute the optimal radius r*(F). Finally, by spanning the 
whole range of interesting fluid flow rates F, we can compute the function relating the mean shear 
rate in the vessel to the vessel radius with the formula: ( 7 )* (r*(F)) = ■ 
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